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We  present  a  study  of  the  collisional  relaxation  of  ion  velocities  in  a  strongly  coupled,  ultracold  neutral 
plasma  on  short  time  scales  compared  to  the  inverse  collision  rate.  The  measured  average  velocity  of  a 
tagged  population  of  ions  is  shown  to  be  equivalent  to  the  ion-velocity  autocon'elation  function.  We  thus 
gain  access  to  fundamental  aspects  of  the  single-particle  dynamics  in  sttongly  coupled  plasmas  and  to  the 
ion  self-diffusion  constant  under  conditions  where  experimental  measurements  have  been  lacking. 
Nonexponential  decay  towards  equilibrium  of  the  average  velocity  heralds  non-Markovian  dynamics 
that  are  not  predicted  by  traditional  descriptions  of  weakly  coupled  plasmas.  This  demonstrates  the 
utility  of  ultracold  neutral  plasmas  for  studying  the  effects  of  strong  coupling  on  collisional  processes, 
which  is  of  interest  for  dense  laboratory  and  astrophysical  plasmas. 
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I.  INTRODUCTION 

In  strongly  coupled  plasmas  [1],  the  Coulomb  interaction 
energy  between  neighboring  particles  exceeds  the  kinetic 
energy,  leading  to  nonbinary  collisions  that  display  tem¬ 
poral  correlations  between  past  and  future  collision  events. 
Such  non-Markovian  dynamics  invalidates  traditional 
theory  for  collision  rates  [2-4]  and  transport  coefficients 
[5,6]  used  for  weakly  coupled  plasmas  and  frustrates  the 
formulation  of  a  tractable  kinetic  theory.  This  challenging 
fundamental  problem  is  also  one  of  the  major  limitations  to 
our  ability  to  accurately  model  equilibration,  transport,  and 
equations  of  state  of  dense  laboratory  and  astrophysical 
plasmas  [7,8],  which  impacts  the  design  of  inertial- 
confinement-fusion  experiments  [9,10],  stellar  chronom- 
etry  based  on  white  dwarf  stars  [11,12],  and  models  of 
planet  formation  [13].  Molecular  dynamics  (MD)  simu¬ 
lations  have  been  the  principal  recourse  for  obtaining  a 
microscopic  understanding  of  short-time  collision  dynam¬ 
ics  in  this  regime  [14-18],  but  direct  comparison  of  results 
with  experiment  has  not  been  possible. 

In  the  experiments  described  here,  we  study  the  effects 
of  strong  coupling  on  collisional  processes  by  measuring 
the  velocity  autocorrelation  function  (VAF)  for  charges  in  a 
strongly  coupled  plasma.  The  VAF,  a  central  quantity  in 
the  statistical  physics  of  many-body  systems,  encodes  the 
influence  of  correlated  collision  dynamics  and  system 
memory  on  individual  particle  trajectories  [19],  and  it  is 
defined  as 
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Z(t)  =  ^(\k(t) -\k( 0)).  (1) 

Here,  \k  is  the  velocity  of  particle  k,  and  brackets  indicate 
an  equilibrium,  canonical-ensemble  average.  Remarkably, 
we  obtain  this  individual-particle  quantity  from  measure¬ 
ment  of  the  bulk  relaxation  of  the  average  velocity  of  a 
tagged  subpopulation  of  particles  in  an  equilibrium  plasma. 
This  contrasts  with  measurements  of  macroscopic-particle 
VAFs  based  on  statistical  sampling  of  individual  trajecto¬ 
ries,  which  is  commonly  used  in  studies  of  dusty-plasma 
kinetics  [20,21]  and  Brownian  motion  [22-24]. 

The  VAF  also  provides  information  on  transport 
processes  since  its  time  integral  yields  the  self-diffusion 
coefficient  through  the  Green-Kubo  relation 


D=  j  Z(t)dt, 

JO 


which  describes  the  long-time  mean-square  displacement 
of  a  given  particle  through  D  =  limf^0O(|r(f)  —  r(0)|2)/6f 
[25].  Our  results  provide  the  first  experimental  measure¬ 
ment  of  the  VAF  and  of  self-diffusion  in  a  three- 
dimensional  strongly  coupled  plasma.  These  results  are 
found  to  be  consistent  with  MD  simulations  to  within  the 
experimental  uncertainty  [5,6]. 

Measurements  are  performed  on  ions  in  an  ultracold 
neutral  plasma  (UCNP),  which  is  formed  by  photoionizing 
a  laser-cooled  atomic  gas  [26,27].  Shortly  after  plasma 
creation,  ions  equilibrate  in  the  strongly  coupled  regime 
with  Coulomb  coupling  parameter, 


T,  = 


4  neQkBTi 


4  trn^j  1  /3 


(2) 
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as  large  as  ~4.  Here,  Tt  is  the  ion  temperature  and  n  is  the 
density.  Electrons  in  the  plasma  provide  a  neutralizing 
background  and  static  screening  on  the  ionic  time  scale 
with  a  Debye  screening  length  A.  This  makes  UCNPs  a 
nearly  ideal  realization  of  a  Yukawa  one-component  plasma 
(OCP)  [28],  a  paradigm  model  of  plasma  and  statistical 
physics  in  which  particles  interact  through  a  pairwise 
Coulomb  potential  screened  by  a  factor  exp(— r/X). 
Ultracold  plasmas  are  quiescent,  near  local  equilibrium, 
and  “clean”  in  the  sense  they  are  composed  of  a  single 
ion  species  and  free  of  strong  background  fields.  Strong 
coupling  is  obtained  at  relatively  low  density,  which 
slows  the  dynamics  and  makes  short-time-scale  processes 
(compared  to  the  inverse  collision  rate)  experimentally 
accessible. 

Powerful  diagnostics  exist  for  dense  laboratory  plasmas. 
However,  their  interpretation  is  complicated  by  the  tran¬ 
sient  and  often  nonequilibrium  nature  of  the  plasmas,  and 
they  do  not  provide  model-independent  information  on  the 
effects  of  particle  correlations  at  short  time  scales.  These 
include,  for  example,  measurement  of  the  dynamic  struc¬ 
ture  factor  with  x-ray  Thomson  scattering  [29-31]  and 
measurement  of  electrical  conductivity  using  a  variety  of 
techniques  [29,32-35].  Comprehensive  studies  of  self- 
diffusivity  [20,21,36-38]  exist  for  strongly  coupled  dusty 
plasmas,  but  these  systems  are  typically  two  dimensional 
and,  therefore,  do  not  directly  illuminate  the  kinetics  of 
bulk,  three-dimensional  plasmas. 


parameter  k  =  a/ A  =  0.1-0.58  in  these  experiments,  for 
Wigner-Seitz  radius  a  =  (3/4/r/;) |  !. 

The  plasma  density  distribution  is  Gaussian  in  shape, 
n  =  »o  exPb  ;'2/2<7(0]’  with  initial  size  (t(0)  =  1-2  mm. 
Because  of  electron  pressure  forces,  the  plasma 
expands  according  to  a 2(f)  =  rr(0)(l  +  f2/Texp)>  where 
rexp  =  10-50  fts  is  the  expansion  time  [44]. 

An  optical  pump-probe  technique  [4,45]  is  used  to 
measure  (A vx(t))+  =  (vx(t))+  -  (vx(t)) totai,  the  average 
velocity  of  a  “spin-tagged”  subpopulation  of  ions  (labeled 
“+”)  relative  to  the  local  bulk  velocity  of  all  the  ions 
((t+(f))total).  Appendix  A  provides  a  proof  that  the  nor¬ 
malized  VAF  'T(f)  =  Z(f)/Z( 0)  is  equivalent  to  the  observ¬ 
able  (Az+(f))+/(Airc(0))+  as  long  as  the  total  system  is 
near  thermal  equilibrium  and  if  terms  beyond  second  order 
are  negligible  in  a  Hermite -Gauss  expansion  of  the  initial 
x- velocity  distribution  function  for  the  subgroup  fx+.  As 
shown  below,  our  experiment  satisfies  these  conditions, 
which  provides  a  new  technique  for  measuring  the  VAF. 

The  evolution  of  (A vx(t))+  is  measured  by  first 
using  optical  pumping  to  create  electron-spin-tagged  ion 
subpopulations  with  nonzero  average  velocity  (Fig.  1). 
Pumping  is  accomplished  by  two  counterpropagating, 
circularly  polarized  laser  beams  each  detuned  by  the  same 
small  amount,  —  Ap,  from  the  2 .S' 1/2  —  2 Pi/ 2  transition  at 
421.7  nm.  Taking  advantage  of  the  unpaired  electron  in  the 
2S 1/2  ground  state,  ions  are  pumped  out  of  the  +1/2  spin 


H.  METHODS 

We  perform  experiments  on  ultracold  neutral  plasmas 
[27],  which  are  created  by  first  laser  cooling  88  Sr  atoms  in  a 
magneto-optical  trap  [39].  Atoms  are  then  photoionized 
with  one  photon  from  a  narrow -band  laser  resonant  with  the 
principal  'So— 'Ri  transition  at  461  nm  and  another 
photon  from  a  tunable  10-ns,  pulsed  dye  laser  near 
413  nm.  The  electron  temperature  ( Te )  in  the  plasma  is 
determined  by  the  excess  photon  energy  above  the  ioniza¬ 
tion  threshold,  which  can  be  tuned  to  set  Te  =  1-1000  K. 
Ions  initially  have  very  little  kinetic  energy,  but  they 
possess  an  excess  of  Coulomb  potential  energy,  and  they 
equilibrate  on  a  microsecond  time  scale  to  a  temperature 
T i  =  0.5-2. 5  K,  determined  primarily  by  the  plasma 
density  [40,41].  The  ion  equilibration  process  is  called 
disorder-induced  heating  (DIH). 

Disorder-induced  heating  limits  the  ions  to  T,  «  2-4.  To 
obtain  measurements  on  more  weakly  coupled  systems 
(r,  <  1),  the  plasma  is  heated  with  ion  acoustic  waves 
[42,43].  Waves  are  excited  by  placing  a  grating  (10  cycles/ 
mm)  in  the  path  of  the  ionization  beam,  which  is  then 
imaged  onto  the  magneto-optical  trap  for  greatest  contrast 
to  create  a  plasma  with  a  striped  density  modulation. 
After  sufficiently  long  time,  the  waves  completely  damp, 
heating  the  plasma  and  reducing  T,.  Electrons  provide  a 
uniform  screening  background  for  the  ions,  with  screening 


FIG.  1.  (a)  Optical  pumping  and  LIF  spectroscopy.  Ions  are 

optically  pumped  from  the  +1/2  to  —1/2  electronic  spin  state 
around  vx  =  —A  p/k,  and  from  the  — 1/2  to +  1/2  spin  state  around 
vx  =  +A p/k  with  two  counterpropagating,  circularly  polarized 
laser  beams  detuned  by  —  Ap  rad/s  from  the  2Sl/2—2Pi/2  tran¬ 
sition  of  the  strontium  ion  (421 .7  nm).  The  velocity  profiles  of  the 
individual  spin  populations  are  measured  with  LIF  using  a  tunable 
circularly  polarized  probe  beam  of  variable  detuning  Apr. 
(b)  Idealized  illustration  of  a  pumped  velocity  distribution  for 
+  1/2  ions  resulting  from  optical  pumping  (red  curve),  along 
with  an  unperturbed  Gaussian  thermal  distribution  (blue  curve). 
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state  and  into  the  —1/2  state  around  the  negative  x  velocity 
vx  =  —Ap/k,  while  ions  are  pumped  from  —1/2  to  +1/2 
spin  around  the  positive  velocity  +A p/k  (the  quantization 
axis  is  taken  to  be  along  the  axis  of  the  pump  beams, 
defined  as  x).  This  creates  subpopulations  of  +1/2  and 
—  1/2  spin  ions  having  velocity  distributions  skewed  in 
opposite  directions,  while  the  entire  plasma  itself  remains 
in  equilibrium.  The  pump  detuning,  Ap/ln  —  30  MHz, 
is  resonant  for  ions  with  \vx\  =  12.6  m/s,  which  is  on 
the  order  of  the  thermal  velocity,  vtb  —  \/kBT ifMi  = 
13.7  m/s  for  Tl  =  2  K. 

We  probe  the  ion  distribution  with  spatially  resolved 
laser-induced  fluorescence  (LIF)  spectroscopy  [44,46] 
(Fig.  1).  A  LIF  probe  beam  tuned  near  the  S 1/2  —  ~P  1/2 
transition  of  the  88Sr  ion  propagates  along  the  x  axis  and 
excites  fluorescence  that  is  imaged  onto  an  intensified  CCD 
camera  with  lx  magnification  (12.5  pm  per  pixel),  from 
which  the  plasma  density  and  x-velocity  distribution  are 
extracted.  By  using  a  circularly  polarized  LIF  probe  beam, 
propagating  nearly  along  the  pump  beam  axis,  we  selec¬ 
tively  probe  only  the  +1/2  ions. 

Pumping  is  applied  several  plasma  periods  (2 n/oop) 
after  ionization  to  allow  the  plasma  to  approach  equili¬ 
brium  after  the  disorder-induced  heating  phase  [27,40]. 
{cop  =  ne2/e0Mj  ~  107  s-1  is  the  ion  plasma  oscillation 
frequency.)  Electro-optic-modulator  pulse  pickers  are  used 
to  achieve  10-ns  time  resolution  for  application  of  the  pump 
and  probe  beams.  The  optical  pumping  time  is  200  ns,  and 
the  pump  intensity  is  200  mW/cm2  (saturation  parameter 
so  =  3.5).  LIF  data  are  taken  at  least  35  ns  after  the  turn-off 
of  the  pump  to  avoid  contamination  of  the  signal  with  light 
from  the  decay  of  atoms  promoted  to  the  2P\/2  state  during 
the  pumping  process.  The  origin  of  time  is  taken  as  the 
middle  of  the  235-ns  preparation  time.  The  pumping  and 
imaging  transition  is  not  closed,  and  about  1/15  of  the 
excitations  result  in  an  ion  decaying  to  a  metastable  2D3,2 
state  that  no  longer  interacts  with  the  lasers.  To  ensure  that 
Larmor  precession  of  the  prepared  atomic  states  does  not 
contaminate  the  data,  a  4.5-G  magnetic  field  is  applied 
along  the  pump-probe  beam  axis. 

The  LIF  spectra  (Fig.  2)  S(v)  are  fit  to  a  convolution 
of  a  Lorentzian  function  L(y)  —  Yl/ {^[(Yl)2  +  2(z/)2]}  of 
frequency  v  (Hz)  with  the  one-dimensional  ion-velocity 
distribution  along  the  laser  axis,  fXt+(vx), 

f  +00 

S(v)  =  /  (1  /X)L(v-s/X)fX'+(s)ds,  (3) 

where  X  is  the  laser  wavelength  and  yl  is  the  width  of 
Lorentzian  spectral  broadening.  For  this  system, 
Yl  =  Yi  +  /,,+  !  +  +),  where  //  is  the  laser  linewidth 
(5.5  MHz)  and  y„  is  the  natural  linewidth  (20.2  MHz). 
The  width  is  power  broadened  by  the  laser  saturation 
parameter  s0-  The  distribution  fx  +  is  modeled  with  a 
Hermite-Gauss  expansion, 


FIG.  2.  (a)-(c)  Time  evolution  of  2Si/2  —  2P\/i  LIF  spectra, 

pumped  (red  circles)  and  unpumped  (green  squares),  for  spin 
+  1/2  ions  in  a  plasma  with  k  =  0.55  and  F,  =  3.0.  (Bars  denote 
time-averaged  values;  see  Appendix  B.)  Spectra  are  fit  (red  and 
green  lines)  to  a  fifth-order  Hermite-Gauss  expansion  of  the 
velocity  distribution  convolved  with  the  Lorentzian  contribution 
from  the  laser  and  natural  linewidths.  (d)-(f)  Full  velocity 
distributions  determined  from  the  fits  of  pumped  data  to  Eqs.  (3) 
and  (4),  along  with  constituent  Hermite  terms  n  =  0-5.  Most  of 
the  velocity  distribution  is  contained  in  the  first  few  terms.  Terms 
of  order  larger  than  zero  decay  in  time  as  the  distribution  relaxes 
to  a  Maxwellian. 
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where  cv  =  \fkBTj Mh  Hn  are  Hermite  polynomials, 
and  the  parameters  C„  are  expansion  coefficients  to  be 
determined  by  a  fit  to  the  data.  For  the  analysis,  N  =  5  is 
chosen  because  the  amplitudes  of  orders  4  and  higher  are 
consistent  with  random  noise. 

The  bulk  velocity  of  the  plasma  (i+(/))total  arises  because 
of  the  plasma  expansion.  To  separate  expansion  from 
velocity  perturbation  due  to  optical  pumping  and  ther¬ 
mal-velocity  spread,  the  plasma  is  divided  into  regions  that 
are  analyzed  independently  [46].  7/  and  (?;A-(r))U)iai  for  each 
region  are  determined  from  analysis  of  LIF  data  from  an 
unpumped  plasma.  LIF  data  from  an  identical  plasma 


021021-3 


T.  S.  STRICKLER  et  al. 


PHYS.  REV.  X  6,  021021  (2016) 


FIG.  3.  Relaxation  of  the  average  velocity  of  spin  +1/2  ions 
in  an  optically  pumped  plasma.  Lines  indicate  memory-function 
and  exponential  fits  to  early-time  data,  as  well  as  a  f-3/2  fit  to  the 
late-time  tail,  (a)  (Ai;x(t))+/(A?;r(0))+  plotted  versus  time  for 
two  strongly  coupled  plasmas,  (b)  Data  from  (a)  versus  time 
scaled  to  the  time  integral  of  aip.  This  plot  shows  the  universal 
scaling  of  (Aiix(t))+/(A?;r(0))+  with  a>p  for  plasmas  of  different 
densities  but  with  approximately  the  same  f ,  and  k.  Panels  (c)  and 
(d)  are  the  same  plots  as  (a)  and  (b)  but  for  a  more  weakly  coupled 
plasma. 

except  with  optical  pumping  are  then  fit  to  Eqs.  (3)  and  (4) 
with  only  the  expansion  coefficients  C„  as  fit  parameters  in 
order  to  determine  (vx(t))+,  the  average  velocity  for  +  ions 
including  expansion  and  the  effects  of  pumping.  From 
the  pumped  and  unpumped  data,  we  find  the  difference 
(A vx{t))+  for  each  of  the  regions,  which  are  averaged 
together  to  arrive  at  final  values,  such  as  in  Figs.  3  and  4. 
A  detailed  example  of  this  analysis  is  provided  in 
Appendix  C. 

Because  of  the  plasma  expansion,  the  density  and  thus 
the  plasma  frequency  a>p  decrease  with  time.  To  account  for 
this,  in  Figs.  3(b),  3(d),  and  4  we  show  the  time  evolution 
of  the  averaged  velocity  as  a  function  of  the  scaled  time 
ts  =  Jq  cop(t')dt',  where  the  evolution  of  mp(t )  is  described 
assuming  self-similar  expansion  of  the  Gaussian  distribu¬ 
tion  [44].  All  quoted  plasma  parameters  indicated  by  a 
symbol  with  a  bar  above  it,  such  as  k  and  f,,  are  time 
averages  over  the  period  of  the  measurement,  as  described 
in  Appendix  B.  The  density  typically  varies  by  a  factor 
of  2  during  the  measurement.  The  normalization  factor 
( A  +(()))  +  is  determined  from  the  fit  of  the  data  using  a 
memory-function  formalism  described  below. 

III.  RESULTS  AND  DISCUSSION 

Sample  LIF  spectra  at  various  times  after  optical  pumping 
are  shown  in  Figs.  2(a)-2(c)  for  a  plasma  with  k  =  0.55 


FIG.  4.  Early-time  behavior  of  (Avx(t)) +/ {Avx(0)) +.  Lines 
indicate  memory-function  and  exponential  fits.  Data  and  fits  for 
plasmas  with  (a)  k  =  0.55  and  f ;  =  3.0  (Te  =  19 K),  (b)  k  =  0.57 
and  f,-  =  3.2  (fe  =  26K),  (c)  k  =  0.24  and  f;  =  2.1 

(Te  =  88 K),  and  (d)  k  =  0.59  and  f ;  =  0.7  ( Te  =  15 K). 
Deviation  from  exponential  decay  is  evident  in  more  strongly 
coupled  plasmas. 


and  f,  =  3.0.  Figures  2(d)-2(f)  show  the  corresponding 
ion-velocity  distributions  and  individual  Hermite-Gauss  com¬ 
ponents  of  the  pumped  velocity  distributions  extracted  from 
fits  to  the  raw  spectra.  At  early  times,  there  is  significant 
amplitude  in  the  n  =  1  term,  corresponding  to  the  skew  in  the 
velocity  distribution.  This  decays  away  as  fx+  approaches  a 
Maxwellian  centered  around  (i+(f))total.  Higher-order  terms 
are  small  at  all  times,  satisfying  an  important  condition  for  the 
proof  that  T»(f)  «  (Avx(t))+/ (Avx(0))+. 

In  the  frame  comoving  with  (t+(f))total  in  a  given  region, 
the  average  velocity  of  +1/2  ions  a  time  t  after  optical 
pumping  is  (A vx(t))+  =  J  fx,+  (vx)(vx  -  ( vx)total)dvx / 
f fx,+{vx)dvx.  Figures  3(a)  and  3(b)  show  sample  data 
for  two  plasmas  of  different  densities  but  with  nearly  the 
same  values  of  k  and  T,  (one  with  k  =  0.55  and  f,  =  3.0, 
the  other  with  k  =  0.57  and  f,  =  3.2).  In  Fig.  3(a), 
data  are  plotted  versus  time,  while  Fig.  3(b)  plots 
(AuJC(f))+/(AuJC(0))+  versus  time  scaled  by  coZ1,  appro¬ 
priately  accounting  for  the  time-changing  value  of  cop. 
The  correspondence  between  the  two  data  sets  apparent  in 
Fig.  3(b)  demonstrates  the  existence  of  a  universal  time 
scale  for  the  dynamics  [28].  Corresponding  data  for 
k  —  0.59,  f ,  =  0.7  are  shown  in  Figs.  3(c)  and  3(d). 

A.  Observation  of  non-Markovian  dynamics 

Data  for  (Avx(t))+/ (Avx(0))+  show  nonexponential 
decay  of  the  average  velocity  up  to  times  given  by 
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f^a/pdf  ~  1,  which  is  a  hallmark  of  non-Markovian 
dynamics  reflecting  the  strong  coupling  of  the  ions.  This 
is  most  clearly  shown  in  Fig.  4,  which  is  an  expanded  view 
of  the  early-time  data  from  velocity  relaxation  curves.  The 
early-time  behavior  of  'P(f)  «  (Avx(t))+/ (A.vx(0))+  can 
be  described  using  a  memory-function  formalism  that  treats 
the  effects  of  collisional  correlations  at  the  microscopic 
level  [15,19,47,48].  It  can  be  derived  from  a  generalized 
Langevin  equation  describing  the  motion  of  a  single  test 
particle  experiencing  memory  effects  and  fluctuating 
forces,  which  is  familiar  from  treatments  of  Brownian 
motion  [19,47].  The  evolution  of  the  VAF  is  found  to  be 


Z(t)  —  -  ['  K(t-  t')Z(t')dt'.  (5) 

Jo 


Here,  K(t  —  t')  is  the  memory  function  describing  the 
influence  at  time  t  from  the  state  of  the  system  at  t' . 

A  general,  closed-form  expression  for  K(t  —  f)  is 
lacking,  but  there  are  expressions  derived  from  simplifying 
assumptions  that  agree  well  with  molecular  dynamics 
simulations  for  simple  fluids  [19,49]  and  Yukawa  potentials 
when  r,  >  20  [14,50].  Some  formulas  introduce  a  time 
constant  r  that  may  be  interpreted  as  the  correlation  time  for 
fluctuating  forces  [15,19,47].  If  one  assumes  that  collisions 
are  isolated  instantaneous  events,  r  ->  0  and  the  memory 
function  becomes  a  delta  function.  This  is  the  Markovian 
limit  in  which  the  evolution  of  a  system  is  entirely 
determined  by  its  present  state  and  T '(f)  has  purely 
exponential  dependence.  Data  from  a  more  weakly  coupled 
sample  [Fig.  4(d),  f ,  =  0.7]  show  no  discernible  rollover 
in  'l'(r)  at  short  times. 

An  often-used  approximation  for  K(t  —  t')  for  the  non- 
Markovian  regime,  valid  for  short  times  and  moderately 
strong  coupling,  is  the  Gaussian  memory  function 
[15,19,47], 


KG(t-f) 


2/c 

'Jim1 


exp 


(t-O2 1 

2t2 


(6) 


which  satisfies  the  conditions  that  memory  effects  vanish 
at  long  time  and  that  it  agrees  with  a  Taylor  expansion 
of  K(t  —  t')  to  second  order  around  t  —  t'  relating  r  to 
frequency  moments  of  the  Fourier  transform  of  the  VAF 
[19].  For  the  Yukawa  OCP,  MD  simulations  have  shown 
that  Eq.  (6)  accurately  reproduces  the  ion  VAF  for  T,  <  10 
and  Wpt  <  n,  and  the  parameter  yc  can  be  related  to  a 
well-defined  collision  rate  [4]. 

Figures  3  and  4  show  fits  of  the  data  in  the  scaled 
time  range  0  <  Jj  copdt'  <  4  to  Eq.  (5)  with  a  Gaussian 
memory  function,  along  with  exponential  fits  of  data  with 
0.8  <  J0r  mpdt'  <  4.5.  At  early  times  ( JJ  a)pdt'  <  0.5),  the 
memory-kernel  fit  captures  the  rollover,  which  is  indicative 
of  non-Markovian  collisional  dynamics.  The  values  of  r 
we  extract  from  the  fit  are  on  the  order  predicted  by  MD 


simulations  of  a  classic  OCP  [15],  although  improved 
experimental  accuracy  is  required  before  a  precise  com¬ 
parison  can  be  made. 


B.  Ion  VAF  and  self-diffusion  coefficients 

Using  (A vx(t))+/ (Avx(0))  +  as  an  approximation  for  the 
normalized  ion  VAF,  the  self-diffusion  coefficient  D  may 
be  calculated  from  our  measurements.  As  is  normally  the 
case  with  calculations  of  this  type,  proper  treatment  of  the 
upper  limit  of  integration  in  the  Green- Kubo  formula  is 
critical  for  obtaining  accurate  results.  The  behavior  of  the 
long-time  tail  of  the  VAF  for  a  Yukawa  system  has  not  been 
explored  in  detail  for  the  regime  of  our  experiment,  and 
this  is  an  important  area  for  future  study.  For  concreteness, 
we  assume  a  f-3/2  dependence,  which  is  well  established 
as  t  ->  co  for  neutral  simple  fluids  [19]  and  is  generally 
accepted  as  the  slowest  possible  decay  [51].  We  fit 
the  last  few  data  points  [beyond  the  time  fcut,  where 
Jo cul  (Op(t')dt'  =  3 to  a  curve,  where  b  is  the  fit 

parameter  and  ts  =  Jj  cop(t')dtl  is  the  scaled  time.  The 
dimensionless  self-diffusion  coefficient,  D*  =  D/a2cop,  is 
thus  calculated  as 


D* 


2b 

Z(ts)dts+^t 


-1/2 

s,N  » 


(7) 


where  the  first  term  is  calculated  numerically  from  the 
data  by  linear  interpolation  and  the  trapezoidal  rule.  The 
time  of  the  last  data  point  is  tN,  and  ts  N  =  JJ  (op{t')dt '  in 
scaled  units.  Extracted  values  for  D*,  along  with  theoretical 
curves  for  k  =  0  and  k  =  0.6  determined  from  a  fit  to 
molecular  dynamics  simulations  [5,6],  are  shown  in  Fig.  5. 
For  comparison.  Fig.  5  also  shows  D*  calculated  from 
the  traditional  Landau-Spitzer  approach  and  from  two 
Chapman-Enskog  calculations.  Of  the  latter,  one  uses  a 
screened  Coulomb  interaction  that  includes  electron  and 
ion  screening  [52,53].  The  other  uses  a  hypernetted  chain 
effective  potential  and  a  modified-Enskog  correction 
calculated  in  Ref.  [54]  (EPT-HNC/Enskog). 

The  error  bars  for  D*  include  statistical  and  systematic 
contributions  we  describe  in  Appendix  B.  The  contribution 
from  the  unmeasured  long-time  tail  of  \D(f)  adds  the 
dominant  uncertainty.  The  lower  error  bar  for  D*  assumes 
no  contribution  beyond  our  measured  points,  while  the 
upper  error  bar  reflects  a  tail  contribution  twice  as  large  as 
our  best  estimate.  This  is  conservative  given  that  the  VAF  is 
exponential  in  the  weakly  coupled  limit.  There  are  signifi¬ 
cant  experimental  improvements  that  can  be  made  in  the 
measurement  and  important  systematic  effects  that  must  be 
investigated.  The  latter  are  scientifically  interesting  in  their 
own  right,  such  as  the  time  scale  for  the  approach  to 
equilibrium  of  velocity  correlations  after  plasma  creation 
and  the  effect  of  plasma  expansion  on  the  microscopic 
dynamics.  Any  complications  caused  by  these  systematics 
can  be  greatly  reduced  by  performing  measurements 
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FIG.  5.  Plot  of  the  normalized  self-diffusion  coefficient  D* 
[Eq.  7)]  calculated  for  our  data.  The  solid  blue  and  solid  red  lines 
represent  results  from  MD  simulations  for  k  =  0  and  k  =  0.6, 
respectively  [5,6].  The  black  dashed  line  represents  the  Landau- 
Spitzer  (LS)  theory  for  weakly  coupled  plasmas,  which  diverges 
in  the  regime  of  strong  coupling.  The  orange  and  the  purple 
dashed  lines  use  screened  Coulomb  and  hypernetted  chain 
interactions,  respectively,  as  described  in  the  text. 

on  a  larger  plasma  for  which  the  expansion  time  scale  is 
much  greater  than  the  characteristic  collisional  time  scale 

(®/Aex p  !)• 

Other  sources  of  uncertainty  include  variation  of  density 
across  the  analysis  regions,  the  spread  in  bulk  plasma 
velocity  across  and  within  regions,  the  time-evolving 
density,  and  the  uncertainty  in  the  density  calibration. 
These  uncertainties  are  reflected  in  the  horizontal  error 
bars  in  Fig.  5  and  can  be  significantly  reduced  in  future 
experiments. 

IV.  CONCLUSIONS 

Utilizing  a  spin-tagging  technique  to  perturb  and  then 
measure  the  average  velocity  (vx(t))+  of  a  subpopulation 
of  ions  in  an  ultracold  neutral  plasma,  and  taking  advantage 
of  an  identification  of  the  time  evolution  of  this  quantity 
with  the  normalized  ion  VAF  [T'(r)],  we  experimentally 
measure  the  VAF  in  a  strongly  coupled  plasma.  From  this, 
we  calculate  the  ion  self-diffusion  coefficient  D,  which 
provides  an  experimental  benchmark  that  has  been  lacking 
for  molecular  dynamics  simulations  of  strongly  coupled 
Coulomb  systems  in  three  dimensions.  The  data  also 
display  a  nonexponential  decay  of  T'(t)  at  early  times, 
which  has  not  previously  been  observed  experimentally  in  a 
bulk  plasma  and  is  indicative  of  non-Markovian  collisional 
dynamics.  This  behavior  is  well  described  by  a  memory- 
function  formalism. 


Overall,  these  measurements  experimentally  validate 
foundational  concepts  describing  how  ion  velocity  corre¬ 
lations  at  the  microscopic  level  determine  the  dynamics 
of  strongly  coupled  systems  at  the  macroscopic  level, 
which  cannot  be  adequately  described  by  simple  analytical 
methods.  Because  ultracold  neutral  plasmas  offer  a  clean 
realization  of  the  commonly  used  Yukawa  OCP  model, 
these  results  are  relevant  for  fundamental  kinetic  theory  and 
other  plasmas  for  which  effects  of  strong  coupling  are 
important. 
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APPENDIX  A:  (Avx  (*))+/ (Avx(0))+  AS 
AN  APPROXIMATION  TO  V{t) 

We  show  that  the  quantity  (At;x(f))+/(Az;x(0))+  we 
measure  in  these  experiments  corresponds  to  the  normal¬ 
ized  VAF  ['k(f)  =  Z(f)/Z(0)]  if  the  initial  velocity  dis¬ 
tribution  for  +1/2  ions  (/+)  is  Maxwellian  in  vy  and  vz  and 
well  described  by  a  second-order  Hermite-Gauss  expansion 
in  vx,  and  the  optical  pumping  prepares  the  subsystem  of 
+  1/2  ions  in  a  nonequilibrium  state  close  to  thermody¬ 
namic  equilibrium. 

To  prove  this,  we  transform  into  the  frame  comoving 
with  any  bulk  hydrodynamic  velocity  of  the  ensemble, 
making  (rr(f))totai  =  0,  which  does  not  invalidate  any  steps 
in  the  proof  but  makes  (A vx(t))+  =  (vx(t))+.  For  simplic¬ 
ity,  we  assume  that  the  plasma  is  spatially  homogeneous 
in  a  volume  V,  although  the  following  arguments  can  be 
readily  extended  to  account  for  nonuniform  distributions. 
Finally,  we  assume  that  the  optical  pumping  occurs 
instantaneously  at  some  initial  time  t  —  0. 

Let  us  consider  a  specified  +1/2  ion  labeled  “s”  with 
position  rs(t)  and  velocity  vs(t)  at  time  t\  for  a  statistical 
description  of  the  dynamics,  it  is  useful  to  define  the 
microscopic  phase  space  density  as  V?(r,  v,  t)  — 
<5(r  —  rs(f))<5(v  —  vs(t))  and  its  statistical  average 
/s.(r,  v,  t).  Before  optical  pumping,  for  t  <  0,  the  system 
is  in  thermal  equilibrium  and  fs(r,\,t)  =  (N(r,  v,  t))  — 
/M(v)/V,  where  /M(v)  is  the  Maxwell-Boltzmann  velocity 
distribution.  After  pumping,  the  subsystem  is  out  of  thermal 
equilibrium: 

/i(r,  v,  t)  =  fM(\)/V  +  <5/s(r,  v,  f).  (Al) 
We  assume  \5fs/{fM/V)\  <SC  1,  so  that  Sfs  satisfies 
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8fs( r,v,  t)=J  dr1  d\'Rs(r  —  r';  v,  v7;  f)5/s(r',  v7, 0), 

(A2) 

where  Rs(r  —  r';  v,  v';  t)  is  the  propagator,  or  retarded 
Green’s  function,  of  the  equation  that  governs  the  temporal 
evolution  of  8fs,  which  is  obtained  by  linearizing  the 
exact  evolution  equation  satisfied  by  fs.  Remarkably, 
the  propagator  Rs,  which  describes  the  nonequilibrium 
dynamics  of  the  system,  is  related  simply  to  the 
equilibrium  time-correlation  function  Cs.(r  —  r';  v,  v';  t)  = 
(Ns( r,  v,  t)Ns( r',  v',  t))  as  follows: 

Cj(r;v,  V;  0  =  /  dr1  J  d\"Rs(r  -  r";  v,  v";  t) 

X  Cs(r';v",  v';0) 

=  Rs(  r;  v,  V;  t)fM(V)/V,  (A3) 

where,  in  the  last  equality,  we  use  the  initial  value 
Cs(r; v,  v';  0)  =  8(r)3(v  —  v')/M(v)/V.  This  relation  is 
an  expression  of  the  fluctuation-dissipation  theorem  [55]. 

From  Eq.  (A2),  the  average  particle  velocity  along  the  x 
direction  determined  in  our  experiment  is 


(»*(*))+  =  J  drd\vx8fs{ r,v,  t) 

=  J  dy  I  dV'V*Rs(y’y,’t)Sfs(y'’0')’ 

where  Rs(\,  v';f)  =  f  drRs(r;v,  V;  t)  and  <5/s(v,  0)  = 
f  drSfs(r,  v,  0).  If,  as  found  experimentally  (see  Fig.  2), 
<5 fs  is  initially  well  described  by  a  second-order  Hermite- 
Gauss  expansion  in  vx,  then 

(»*(*)>+  =  J  dy  J  d\'vx[c0  +  c1v'x  +  c2(v,x)2] 

xRs(v,  y'-,t)fM(V)/V 

=  j  dy  f  dy’Vxic0  +  cl  vx  +  C2(v'x)2} 

X  Cs(v,  v';  /), 

where  we  use  the  fluctuation-dissipation  theorem  Eq.  (A3) 
to  obtain  the  second  equality.  The  zeroth-  and  second-order 
terms  vanish  since  Cs(v,  v';  f)  =  Cs(— v,  — v';  f),  and  the 
previous  result  simplifies  to 

(vx(t))+  =  c i  J  d\  J  d\'vxv'xCs{y,\'\t) 

=  Ci(vXjS(t)vX'S(0))  =  c{Z(t). 

Therefore,  we  obtain  the  desired  result: 


(g(0)+ 

<®*(0))+ 


Z(t) 

Z(0) 


(A4) 


APPENDIX  B:  VARIATION  AND 
UNCERTAINTIES  OF  PLASMA  PARAMETERS 

Because  the  plasmas  in  these  experiments  undergo 
expansion  during  the  measurements,  the  values  of  various 
plasma  parameters  quoted  for  a  given  VAF  measurement 
are  time-averaged  values  given  by  a  trapezoidal  approxi¬ 
mation  to  @  =  [  f()J  Q(t')dt']/tf,  where  tr  is  the  total  time  of 
measurement  and  ©(f)  represents  T,(f),  n,(f),  /c(f),  cop(t), 
or  r,(f). 

We  now  discuss  the  effects  of  expansion  on  the  values  of 
the  various  plasma  parameters,  as  well  as  the  methods  used 
to  incorporate  these  systematic  variations  into  the  quoted 
uncertainties.  UCNPs  are  created  by  ionizing  neutral  atoms 
in  a  cloud  with  a  3D  Gaussian  profile.  Immediately  after 
ionization,  T,  equals  the  temperature  of  the  trapped  atoms 
(around  10  mK),  but  on  the  time  scale  of  the  inverse  plasma 
frequency,  the  ions  heat  to  a  few  degrees  kelvin  due  to 
disorder-induced  heating.  Measurements  of  the  velocity 
autocorrelation  function  and  the  self-diffusion  constant  are 
performed  after  the  DIH  is  complete  and  the  plasma  is  near 
local  thermal  equilibrium.  After  DIH,  the  ion  temperature 
still  changes  due  to  electron-ion  heating  and  cooling  from 
expansion  into  the  vacuum  [56],  but  this  is  a  slow  evolution 
compared  to  the  experimental  time  scale. 

After  ionization,  the  plasma  expands  due  to  electron 
pressure  [44].  The  expansion  of  the  plasma  cloud  is  self¬ 
similar  and  well  described  by  o(t)2  =  (t(0)2(1  +  /2/reXp), 
where  o(t)  is  the  Gaussian  width  parameter  at  time  t  after 
ionization  and  rexp  is  the  expansion  time  scale  defined  by 

rexp  =  y/mj(y(0)2/kB[Te(0)  +  T,(0)],  where  m,-  is  the  ion 
mass,  kg  is  Boltzmann’s  constant,  7^(0)  is  the  initial 
electron  temperature,  and  Tt( 0)  is  the  ion  temperature 
after  DIH.  Because  of  this  self-similar  expansion,  an 
infinitesimal  spatial  region  of  the  plasma  that  starts 
at  a  given  normalized  radial  position  f (0)  =  rj  remains 
centered  at  r  = ;/  throughout  the  expansion,  where 
?(f)  =  r {t)/a{t).  To  track  a  well-defined  population  of 
ions  in  the  plasma,  the  analysis  regions  move  with  the 
plasma  expansion  at  a  fixed  r ,  while  also  increasing  in  size 
at  the  same  rate  as  a{t). 

The  density  is  determined  from  the  amplitude  of  the  LIF 
spectrum  [46],  with  initial  calibration  given  by  the  DIH 
curve  as  described  in  Ref.  [28].  The  density  at  a  position 
of  constant  f  changes  by  as  much  as  a  factor  of  2  during 
the  course  of  a  measurement,  but  cop  and  T,  vary  as  « 1/2  and 
n 1  ,/3 ,  respectively,  leading  to  a  much  smaller  change  in 
these  quantities. 

The  gradient  of  velocity  outward  from  the  center  of  the 
plasma  broadens  the  LIF  spectrum  in  proportion  to  the 
length  of  the  analysis  region  along  the  direction  of 
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FIG.  6.  Typical  evolution  of  (a)  Th  (b)  nh  (c)  a>p,  (d)  r;,  (e)  k,  and  (f)  Te  versus  scaled  time  with  smooth  fits. 


propagation  of  the  LIF  laser,  which  complicates  the 
determination  of  plasma  temperature.  Small  regions  are 
used  (corresponding  to  50  pi m  along  the  k  vector  of  the 
LIF  beam),  so  that  the  velocity  change  over  a  region  is  less 
than  10%  of  the  thermal  velocity  at  the  very  latest  times. 
During  most  of  the  measurement,  it  is  significantly  less. 

The  effects  of  expansion  are  illustrated  for  a  plasma  with 
k  =  0.24  and  f ,  =  2.1  in  Fig.  6.  This  plasma  has  one  of  the 
fastest  expansions  in  the  data  set,  and  represents  a  “worst 
case”  for  the  complications  introduced  by  plasma  expan¬ 
sion.  In  Fig.  6,  the  averages  of  T,  and  nt  for  all  the  analysis 
regions  are  shown  with  vertical  error  bars  representing  the 
standard  deviation.  cop  and  T,  are  calculated  from  7’,  and  nt 
for  each  region  and  then  averaged.  Data  are  fit  to  smooth 
functions  in  time,  which  are  used  to  calculate  a>p(t)  in  the 
expressions  for  scaled  times  in  Figs.  2  and  3,  and  to 
determine  the  values  of  time-averaged  f ,  used  to  plot  the 
self-diffusion  coefficient  data  in  Fig.  5.  The  smooth  lit  to  T, 
is  a  second-order  polynomial,  while  the  fits  to  the  remain¬ 
ing  quantities  are  made  using  the  plasma  expansion 
formulas  [44].  The  initial  electron  temperature  77(0)  is 
determined  by  the  detuning  of  the  photoionization  laser 
above  threshold  [27],  and  the  evolution  77(f)  is  assumed  to 
follow  predictions  for  a  self-similar  plasma  expansion  [44]. 

The  spread  in  the  value  of  T,  during  a  single  measure¬ 
ment  of  the  self-diffusion  constant,  as  shown  in  Fig.  6,  is 
taken  as  the  systematic  uncertainty  of  f ,  in  Fig.  5,  and  is 
reflected  in  the  error  bars.  Statistical  uncertainties  are  also 


reflected  in  the  error  bars,  but  they  are  a  minor  contribution. 
The  spread  in  mp  values  is  also  taken  as  the  systematic 
uncertainty  of  scaled  time  ts,  which  propagates  into 
the  calculation  of  the  self-diffusion  constant  through  the 
Green-Kubo  integral  of  ( Avx(t))+ / (A?;;t(0))+  over  the 
scaled  time.  Overall,  the  estimate  of  the  total  uncertainties 
in  f ,  and  D  is  a  conservative  one,  incorporating  statistical 
uncertainty,  systematic  uncertainty  arising  from  incomplete 
knowledge  of  the  form  of  the  long-time  tail  of  the  VAF,  and 
“uncertainty”  introduced  by  the  range  of  values  the  plasma 
parameters  take  over  the  course  of  the  measurement. 

APPENDIX  C:  SPECTRA  FOR  INDIVIDUAL 
REGIONS  AND  VARIATION  ACROSS 
THE  PLASMA 

In  Fig.  7(a),  spectra  from  individual  analysis  regions  are 
shown  for  the  plasma  described  in  Fig.  6.  The  unpumped 
spectra  exhibit  a  shift  of  the  center  frequency  from  zero  that 
increases  with  time  and  distance  from  the  plasma  center 
(region  0),  reflecting  the  plasma  expansion.  By  fitting  the 
unpumped  spectra  to  a  Gaussian  distribution  convolved 
with  the  laser  and  natural  linewidth  broadening  (assumed 
to  be  a  Lorentzian  function),  the  center  frequency  and 
Gaussian  width  are  extracted.  The  center  frequency  deter¬ 
mines  the  bulk  expansion  velocity  vx  for  each  region.  The 
width  determines  the  ion  temperature  T T ,  and  vx  are 
set  as  constants  in  fits  of  the  corresponding  pumped 
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FIG.  7.  (a)  Pumped  and  unpumped  LIF  spectra  from  regions  of  width  0.1  cx  across  the  plasma  cloud  spanning  x  =  ±0.2<ta.  for  the 

plasma  described  in  Fig.  6.  (b)  Pumped  and  unpumped  velocity  distributions  determined  from  the  spectra  shown  in  (a). 


distributions  to  the  fifth-order  Hermite  expansion  [Eq.  (5)] 
convolved  with  the  Lorentzian  broadening.  As  Fig.  7(a) 
shows,  the  fit  to  the  raw  data  captures  the  “skewing”  of 
pumped  distributions  well.  The  convolution  of  the  distri¬ 
bution  with  the  Lorentzian  broadening  preserves  the  first 
moment  of  the  distributions,  so  this  broadening  should  not 
affect  the  value  of  the  average  velocity. 

In  Fig.  7(b),  the  underlying  velocity  distributions 
from  the  fits  to  each  raw  spectrum  in  Fig.  7(a)  are  shown, 
along  with  each  term  of  the  Hermite  expansion.  Together, 
Figs.  7(a)  and  7(b)  demonstrate  that  the  first-order  Hermite 
term  dominates  over  the  higher-order  odd  terms  at  early 
times  across  all  the  regions,  which  is  a  requirement  for 


(Avx(t))+/ (Avx(0))+  to  be  identified  as  the  normalized 
velocity  autocorrelation  function. 

To  disentangle  the  effects  of  plasma  expansion  from 
thermal  ion  motion  in  the  LIF  spectra,  it  is  important  to 
use  regions  that  are  relatively  small  along  the  direction  of 
propagation  of  the  LIF  laser  beam.  Seven  small,  overlapping 
regions  are  analyzed  individually.  The  regions  have  widths 
0.1  a  along  the  direction  of  propagation  of  the  LIF  laser  and 
span  from  —0.2ax  to  +0.2ffx  along  the  same  direction. 

The  values  of  (A  vx(t))+  for  each  of  the  seven  regions  are 
averaged  together  to  get  one  time  evolution  of  (A vx(t))  + 
for  the  plasma.  A  typical  result  is  shown  in  Fig.  8  for  the 
plasma  described  in  Fig.  6.  Black  points  are  the  average  of 
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FIG.  8.  Evolution  of  average  velocity  of  pumped  spectra  across 
all  regions  for  the  plasma  described  in  Fig.  6.  The  black  circles 
represent  the  mean  value  for  each  time  point.  The  inset  represents 
a  closeup  view  of  the  same  data  over  a  scaled  time  ranging 
from  0  to  4.5  to  highlight  the  nonexponential  decay  at  early  times. 


all  regions  at  each  time.  The  error  bars  are  the  statistical, 
one-standard-deviation  uncertainties  of  the  mean  values, 
which  we  use  in  the  main  article  to  calculate  the  uncer¬ 
tainties  in  (Avx(t))+/ (Avx(0))+.  The  fits  to  the  exponen¬ 
tial  function  (using  data  in  the  range  0.8  <  copt  <  4)  and  to 
the  memory  kernel  (using  data  in  the  range  0  <  a>pt  <  4) 
are  also  shown  for  reference. 
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